1) Contextualização

O Brasil destaca-se como o maior exportador de café no mundo e o segundo maior consumidor dessa bebida. Responsável por um terço da produção global, o Brasil mantém a liderança na produção de café há mais de 150 anos.

Entre as regiões produtoras, Minas Gerais é o maior estado produtor, contribuindo com cerca de 50% da produção nacional. Já São Paulo - conhecido pela sua tradição no cultivo de café, - produz exclusivamente Arábica.

Por ser uma paixão nacional, também somos o segundo maior consumidor no mundo, perdendo apenas para os EUA.

Em 2023 pela primeira vez na história, um produtor brasileiro, a São Mateus Agropecuária, conquistou o título “Best of the Best” no Prêmio Internacional Ernesto Illy. Esse reconhecimento é resultado de décadas de investimentos em técnicas sustentáveis e agricultura regenerativa. A fazenda foi pioneira em práticas como reuso da água e cuidado com o ecossistema, contribuindo para a qualidade do café brasileiro.

Fonte1: CONAB. Disponível em: https://www.abic.com.br/tudo-de-cafe/o-cafe-brasileiro-na-atualidade/. Acesso em: 02/12/23. Fonte2: Ansa. Disponível em: https://www.abic.com.br/tudo-de-cafe/o-cafe-brasileiro-na-atualidade/. Acesso em: 02/12/23. Fonte3: Disponível em: https://satocomunicacao.com.br/cafe-e-paixao-nacional-e-segunda-bebida-mais-consumida-pelos-brasileiros/. Acesso em: 02/12/23.

2) Objetivo

O objetivo principal deste projeto é desenvolver e validar um modelo de CLASSIFICÃO, capaz de distinguir cafés brasileiros de não-brasileiros com base em avaliações detalhadas fornecidas por Q Graders profissionais.

3) Início Projeto

3.1) Carregando os pacotes

library(tidymodels)
library(tidyverse)
library(skimr)
library(ggplot2)
library(gridExtra)
library(patchwork)
library(reshape2)
library(car)
library(corrplot)
library(ranger)
library(MASS)
library(keras)
library(rsample)
library(yardstick)
library(FactoMineR)
library(factoextra)
library(doParallel)
library(caret)
library(pROC)
library(purrr)
library(gt)
library(ggrepel)
library(dplyr)
library(rmarkdown)

3.2) Importe das bases

#2)Import da base
dados_originais <- read.csv("coffee_ratings_1.csv")

4) Análise Exploratória

4.1) Conhecendo a base Coffe_Ratings

paged_table(dados_originais)

4.2) Estrutura

# Mostre o número de linhas e colunas do DataFrame
cat("Número de linhas:", nrow(dados_originais), "\n")

## Número de linhas: 1339

cat("Número de colunas:", ncol(dados_originais), "\n")

## Número de colunas: 43

4.3) Contagem de Vazios

#Verificando os dados missing
resultado <- dados_originais %>%
  lapply(type_sum) %>%
  as_tibble() %>%
  pivot_longer(cols = 1:ncol(dados_originais),
               names_to = "Coluna",
               values_to = "Tipo") %>%
  mutate(Contagem_NA = colSums(is.na(dados_originais))) %>%
  arrange(desc(Contagem_NA))

resultado_df <- as.data.frame(resultado)

print(resultado_df)

##                   Coluna Tipo Contagem_NA
## 1             lot_number  chr        1063
## 2              farm_name  chr         359
## 3                   mill  chr         315
## 4               producer  chr         231
## 5    altitude_low_meters  dbl         230
## 6   altitude_high_meters  dbl         230
## 7   altitude_mean_meters  dbl         230
## 8               altitude  chr         226
## 9                variety  chr         226
## 10                 color  chr         218
## 11               company  chr         209
## 12     processing_method  chr         170
## 13            ico_number  chr         151
## 14                region  chr          59
## 15          harvest_year  chr          47
## 16                 owner  chr           7
## 17               owner_1  chr           7
## 18     country_of_origin  chr           1
## 19               quakers  int           1
## 20      total_cup_points  dbl           0
## 21               species  chr           0
## 22        number_of_bags  int           0
## 23            bag_weight  chr           0
## 24    in_country_partner  chr           0
## 25          grading_date  chr           0
## 26                 aroma  dbl           0
## 27                flavor  dbl           0
## 28            aftertaste  dbl           0
## 29               acidity  dbl           0
## 30                  body  dbl           0
## 31               balance  dbl           0
## 32            uniformity  dbl           0
## 33             clean_cup  dbl           0
## 34             sweetness  dbl           0
## 35         cupper_points  dbl           0
## 36              moisture  dbl           0
## 37  category_one_defects  int           0
## 38  category_two_defects  int           0
## 39            expiration  chr           0
## 40    certification_body  chr           0
## 41 certification_address  chr           0
## 42 certification_contact  chr           0
## 43   unit_of_measurement  chr           0

4.4) Dropando Colunas

# Dropando preditoras desinteressantes apos definiciao grupo
dados_tratados_1 <- dados_originais %>% 
  dplyr::select(- moisture,-cupper_points ,-category_two_defects ,-quakers ,-category_one_defects,-owner, -region, -farm_name, -lot_number, -mill, -ico_number, -company, -altitude, -producer, -number_of_bags, -bag_weight,
                -in_country_partner, -harvest_year, -grading_date, -owner_1, -variety, -processing_method, -color,
                -expiration, -certification_body, -certification_address, -certification_contact,-unit_of_measurement,-altitude_low_meters, -altitude_high_meters,-altitude_mean_meters)

4.5) Dropando Dados Inconsistêntes

#4.1) Dropando linhas missing da coluna region
dados_tratados_2 <- dados_tratados_1 %>% 
  dplyr::filter(country_of_origin != "Cote d?Ivoire")

# Excluindo o outlier em Honduras
dados_tratados_4 <- dados_tratados_2 %>%
  filter(!(country_of_origin == "Honduras" & total_cup_points == 0))

5) Análise Descritiva

5.1) Análise Estatistica + Checagem das Frequências

#Estatistica
glimpse(dados_tratados_4)

## Rows: 1,336
## Columns: 12
## $ total_cup_points  <dbl> 90.58, 89.92, 89.75, 89.00, 88.83, 88.83, 88.75, 88.…
## $ species           <chr> "Arabica", "Arabica", "Arabica", "Arabica", "Arabica…
## $ country_of_origin <chr> "Ethiopia", "Ethiopia", "Guatemala", "Ethiopia", "Et…
## $ aroma             <dbl> 8.67, 8.75, 8.42, 8.17, 8.25, 8.58, 8.42, 8.25, 8.67…
## $ flavor            <dbl> 8.83, 8.67, 8.50, 8.58, 8.50, 8.42, 8.50, 8.33, 8.67…
## $ aftertaste        <dbl> 8.67, 8.50, 8.42, 8.42, 8.25, 8.42, 8.33, 8.50, 8.58…
## $ acidity           <dbl> 8.75, 8.58, 8.42, 8.42, 8.50, 8.50, 8.50, 8.42, 8.42…
## $ body              <dbl> 8.50, 8.42, 8.33, 8.50, 8.42, 8.25, 8.25, 8.33, 8.33…
## $ balance           <dbl> 8.42, 8.42, 8.42, 8.25, 8.33, 8.33, 8.25, 8.50, 8.42…
## $ uniformity        <dbl> 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.…
## $ clean_cup         <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, …
## $ sweetness         <dbl> 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 9.3…

skim(dados_tratados_4)

Data summary
Name	dados_tratados_4
Number of rows	1336
Number of columns	12
_______________________
Column type frequency:
character	2
numeric	10
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
species	0	1	7	7	0	2	0
country_of_origin	0	1	4	28	0	35	0

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
total_cup_points	1	82.16	2.69	59.83	81.17	82.50	83.67	90.58	▁▁▁▇▁
aroma	1	7.57	0.32	5.08	7.42	7.58	7.75	8.75	▁▁▂▇▁
flavor	1	7.53	0.34	6.08	7.33	7.58	7.75	8.83	▁▂▇▃▁
aftertaste	1	7.41	0.35	6.17	7.25	7.42	7.60	8.67	▁▃▇▂▁
acidity	1	7.54	0.32	5.25	7.33	7.58	7.75	8.75	▁▁▃▇▁
body	1	7.52	0.31	5.08	7.33	7.50	7.67	8.58	▁▁▁▇▁
balance	1	7.52	0.35	5.25	7.33	7.50	7.75	8.75	▁▁▃▇▁
uniformity	1	9.84	0.49	6.00	10.00	10.00	10.00	10.00	▁▁▁▁▇
clean_cup	1	9.84	0.72	0.00	10.00	10.00	10.00	10.00	▁▁▁▁▇
sweetness	1	9.86	0.55	1.33	10.00	10.00	10.00	10.00	▁▁▁▁▇

#frequencia
table(dados_tratados_4$country_of_origin)

## 
##                       Brazil                      Burundi 
##                          132                            2 
##                        China                     Colombia 
##                           16                          183 
##                   Costa Rica                      Ecuador 
##                           51                            3 
##                  El Salvador                     Ethiopia 
##                           21                           44 
##                    Guatemala                        Haiti 
##                          181                            6 
##                     Honduras                        India 
##                           52                           14 
##                    Indonesia                        Japan 
##                           20                            1 
##                        Kenya                         Laos 
##                           25                            3 
##                       Malawi                    Mauritius 
##                           11                            1 
##                       Mexico                      Myanmar 
##                          236                            8 
##                    Nicaragua                       Panama 
##                           26                            4 
##             Papua New Guinea                         Peru 
##                            1                           10 
##                  Philippines                       Rwanda 
##                            5                            1 
##                       Taiwan Tanzania, United Republic Of 
##                           75                           40 
##                     Thailand                       Uganda 
##                           32                           36 
##                United States       United States (Hawaii) 
##                           10                           73 
##  United States (Puerto Rico)                      Vietnam 
##                            4                            8 
##                       Zambia 
##                            1

table(dados_tratados_4$species)

## 
## Arabica Robusta 
##    1308      28

5.2) Dropando Desbalanceamento

#Grande desbalanceamento de dados entre robusta e Arabica causando efeito de multicolinearidade
dados_tratados_4 <- dplyr::select(dados_tratados_4, -species)

6) Análise Grafica

6.1) Contagem por Pais e ponderação %

# Balanceamento paises

# Tiblbe com a frequencia dos dados Paises 
tabela_regioes <- dados_tratados_4 %>%
  count(country_of_origin) %>%
  rename(Contagem = n) %>%
  group_by(country_of_origin) %>%
  mutate(Percentual_Participacao = (Contagem / nrow(dados_tratados_4)) * 100) %>%
  ungroup() %>%
  arrange(desc(Contagem))

# Adicionar uma linha para somar os totais de contagem e participação
tabela_regioes <- tabela_regioes %>%
  bind_rows(
    tibble(
      country_of_origin = "Total",
      Contagem = sum(tabela_regioes$Contagem),
      Percentual_Participacao = sum(tabela_regioes$Percentual_Participacao)
    )
  )

# Removendo a linha do 'Total' para o gráfico (se não quiser incluir o total)
tabela_regioes <- tabela_regioes %>% filter(country_of_origin != "Total")

# Reordenando os países com base na Contagem para exibição decrescente
tabela_regioes <- tabela_regioes %>% 
  mutate(country_of_origin = reorder(country_of_origin, -Contagem))

# Calculando o máximo da Contagem e Percentual_Participacao para a escala
max_contagem <- max(tabela_regioes$Contagem)
max_percentual <- max(tabela_regioes$Percentual_Participacao)

# Criando o gráfico
grafico <- ggplot(tabela_regioes, aes(x = country_of_origin, y = Contagem)) +
  geom_bar(stat = "identity", fill = "blue") +
  geom_line(aes(y = Percentual_Participacao * max_contagem / max_percentual, group = 1), 
            color = "red", linewidth = 0.5) +
  geom_point(aes(y = Percentual_Participacao * max_contagem / max_percentual), 
             color = "red", size = 2) +
  geom_text(aes(label = Contagem, y = Contagem), vjust = -0.5, size = 3.3) +
  geom_text(aes(y = Percentual_Participacao * max_contagem / max_percentual, label = ifelse(Percentual_Participacao > 1, paste0(round(Percentual_Participacao, 0), "%"), "")), 
            color = "white", vjust = 2.5, size = 2.3) +
  scale_y_continuous(sec.axis = sec_axis(~ . * max_percentual / max_contagem, name = "Percentual de Participação (%)")) +
  labs(x = "Países", y = "Contagem") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Exibindo o gráfico
print(grafico)

#### 6.2) Nota Média por País

# 2.2) Nota media por pais

# Calculando a média de total_cup_points  por país
media_pontuacao_por_pais <- dados_tratados_4 %>%
  group_by(country_of_origin) %>%
  summarise(media_pontuacao = mean(total_cup_points , na.rm = TRUE)) %>%
  arrange(desc(media_pontuacao))

# Criando o gráfico com ggplot2
plot <- ggplot(media_pontuacao_por_pais, aes(x = reorder(country_of_origin, -media_pontuacao), y = media_pontuacao)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = sprintf("%.2f", round(media_pontuacao, 2))),
            position = position_dodge(width = 0.9), hjust = 1.1, size = 3.5, color = "white") +
  theme_minimal() +
  labs(title = "Média de Pontuação Total do Café por País", x = "País", y = "Média de Pontuação") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5),
        plot.title = element_text(size = 18, face = "bold"),
        axis.title = element_text(size = 16)) +
  coord_flip() # Usado para inverter as coordenadas e facilitar a leitura

# Plotando o gráfico
print(plot)

6.3) Dipersão Total Notas vs Pais

# Criar um único gráfico de dispersão para Total Cup Points vs paises 
g <- ggplot(dados_tratados_4, aes(x = country_of_origin, y = total_cup_points)) + 
  geom_point(position = position_jitter(width = 0.2), alpha = 0.5) + 
  geom_smooth(aes(group = country_of_origin), method = "lm", se = FALSE, color = "blue") +
  labs(title = "Dispersao do Total Cup Points por País") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Exibindo o gráfico
print(g)

## `geom_smooth()` using formula = 'y ~ x'

6.4) Comparação das Notas do Brasil vs Países

por Atributo

# 2.2) Nota media Brasil por atributo

# Filtrar os dados para incluir apenas as entradas do Brasil
dados_graf_brasil <- dados_tratados_4 %>%
  filter(country_of_origin == "Brazil") 

# Calcular a média de cada atributo para o Brasil
medias_brasil <- dados_graf_brasil %>%
  summarise(
    aroma = mean(aroma, na.rm = TRUE),
    flavor = mean(flavor, na.rm = TRUE),
    aftertaste = mean(aftertaste, na.rm = TRUE),
    acidity = mean(acidity, na.rm = TRUE),
    body = mean(body, na.rm = TRUE),
    balance = mean(balance, na.rm = TRUE),
    uniformity = mean(uniformity, na.rm = TRUE),
    clean_cup = mean(clean_cup, na.rm = TRUE),
    sweetness = mean(sweetness, na.rm = TRUE)
  ) %>%
  mutate(country = "Brasil")

# Filtrar os dados para incluir apenas as entradas do Brasil
dados_graf_outros <- dados_tratados_4 %>%
  filter(country_of_origin != "Brazil") 

# Calcular a média de cada atributo para os outros países
medias_outros_paises <- dados_graf_outros %>%
  summarise(
    aroma = mean(aroma, na.rm = TRUE),
    flavor = mean(flavor, na.rm = TRUE),
    aftertaste = mean(aftertaste, na.rm = TRUE),
    acidity = mean(acidity, na.rm = TRUE),
    body = mean(body, na.rm = TRUE),
    balance = mean(balance, na.rm = TRUE),
    uniformity = mean(uniformity, na.rm = TRUE),
    clean_cup = mean(clean_cup, na.rm = TRUE),
    sweetness = mean(sweetness, na.rm = TRUE)
  ) %>%
  mutate(country = "Outros Países")

# Unir os dados do Brasil com os dados dos outros países
dados_comparativos <- rbind(medias_brasil, medias_outros_paises) %>%
  gather(attribute, mean_value, -country)

# Criando o gráfico com ggplot2
plot_comparativo <- ggplot(dados_comparativos, aes(x = attribute, y = mean_value, fill = country)) +
  geom_col(position = position_dodge(), width = 0.8) +
  geom_text(aes(label = sprintf("%.2f%%", mean_value)),
            position = position_dodge(width = 0.8), vjust = 1.5, size = 3.7, angle = 45) +
  theme_minimal() +
  labs(title = "Comparação Brasil vs. Outros Países", 
       x = "Atributo", y = "Média") +
  scale_fill_manual(values = c("Brasil" = "steelblue", "Outros Países" = "grey")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        plot.title = element_text(size = 18, face = "bold"),
        axis.title = element_text(size = 16))

# Exibir o gráfico
print(plot_comparativo)

#### 6.5) Dipersão Total Notas vs Pais

# Criar um único gráfico de dispersão para Total Cup Points vs paises 
g <- ggplot(dados_tratados_4, aes(x = country_of_origin, y = total_cup_points)) + 
  geom_point(position = position_jitter(width = 0.2), alpha = 0.5) + 
  geom_smooth(aes(group = country_of_origin), method = "lm", se = FALSE, color = "blue") +
  labs(title = "Dispersao do Total Cup Points por País") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Exibindo o gráfico
print(g)

## `geom_smooth()` using formula = 'y ~ x'

7) Ajuste Final para a Modelagem

7.1) Drop Final ajuste

dados_tratados_4 <- dplyr::select(dados_tratados_4, -total_cup_points)

7.2) Comparação da Base Antes de Depois

#5) Checando as limpezas de dados antes e depois

dim_antes <- dim(dados_originais)
dim_depois <- dim(dados_tratados_4)

# Criando uma tibble para exibir as dimensões
dimensoes <- tibble(
  "Checagem" = c("Antes", "Depois"),
  "Linhas" = c(dim_antes[1], dim_depois[1]),
  "Colunas" = c(dim_antes[2], dim_depois[2])
)
# Exibindo a tibble
print(dimensoes)

## # A tibble: 2 × 3
##   Checagem Linhas Colunas
##   <chr>     <int>   <int>
## 1 Antes      1339      43
## 2 Depois     1336      10

7.3) Comparação da Base Antes de Depois

dados_final <- dados_tratados_4
# Exibindo a tibble
paged_table(dados_final)

8) Preparação para Modelagem

# Definindo a semente para reprodutibilidade
set.seed(15)

Aqui precisamos usar uma técnica de Estratificacao para poder balancear os dados , ja que o número de observações para o Brasil estao desproporcionais comparado com os demais paises. Sem esse balanceamento o modelo ficará tendencioso, ou seja , ele será capaz de classificar com muita precisào o que nào e café brasileiro do que brasileiro , nao sendo esse nosso objetivo.

8.1) Estratificação

# Subamostragem dos dados de outros países para equilibrar com os dados do Brasil
dados_brazil <- dados_final %>% filter(country_of_origin == "Brazil")
dados_outros <- dados_final %>% filter(country_of_origin != "Brazil")

n_brazil <- nrow(dados_brazil)
dados_outros_subamostrados <- dados_outros %>%
  slice_sample(n = n_brazil)

# Combinando os dados subamostrados com os dados do Brasil
dados_combinados <- bind_rows(dados_brazil, dados_outros_subamostrados)

8.2) Treino e Teste

# Redefinindo a divisão de treino e teste com os dados balanceados
split_combinado <- initial_split(dados_combinados, prop = 0.7)
treinamento <- training(split_combinado)
teste <- testing(split_combinado)

8.3) Checagem do Balanceamento - Estratificacao

# Checando o Balanceamento____________________________________________________
n_brasil_treinamento <- treinamento %>%
  filter(country_of_origin == "Brazil") %>%
  summarize(count = n()) %>%
  pull(count)

# Contando os registros do Brasil no conjunto de teste
n_brasil_teste <- teste %>%
  filter(country_of_origin == "Brazil") %>%
  summarize(count = n()) %>%
  pull(count)

# Criando uma tibble para mostrar os resultados
resumo_brasil <- tibble(
  Conjunto = c("Treinamento", "Teste"),
  Total_Brasil = c(n_brasil_treinamento, n_brasil_teste)
)
# Exibindo o resumo final
print(resumo_brasil)

## # A tibble: 2 × 2
##   Conjunto    Total_Brasil
##   <chr>              <int>
## 1 Treinamento           95
## 2 Teste                 37

print(nrow(treinamento))

## [1] 184

print(nrow(teste))

## [1] 80

8.4) Transformando a Coluna Pais em Binária para Classificacao

# Checando o Balanceamento____________________________________________________
# Criar a nova variável alvo para classificação binária como fator
treinamento$from_brazil <- as.factor(ifelse(treinamento$country_of_origin == "Brazil", 1, 0))
teste$from_brazil <- as.factor(ifelse(teste$country_of_origin == "Brazil", 1, 0))

# Remover a coluna original 'country_of_origin'
treinamento$country_of_origin <- NULL
teste$country_of_origin <- NULL

8.5) Pré Processamento - Receita , Prep e Cake

# Preparando a receita com a nova variável alvo
receita <- recipe(from_brazil ~ ., data = treinamento) %>% 
  step_normalize(all_numeric())

receita_prep <- prep(receita)
treinamento_proc <- bake(receita_prep, new_data = NULL)
teste_proc <- bake(receita_prep, new_data = teste)

9) Modelagem

9.1) Logistica

fit_glm <- logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")
logistica_simples <- fit(fit_glm, from_brazil ~ ., data = treinamento_proc)

# Fazendo previsões no conjunto de teste e combinando com os dados de teste
logistica <- predict(logistica_simples, teste_proc, type = "class") %>%
  bind_cols(teste_proc) %>%
  rename(Predito = .pred_class)

# Verificando se a coluna 'from_brazil' está presente
if (!"from_brazil" %in% names(logistica)) {
  logistica$from_brazil <- teste$from_brazil
}

# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
logistica <- logistica %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))

9.1) Ridge Tunado

 #Definindo o modelo de regressão logística com penalidade Ridge
ridge <- logistic_reg(penalty = tune(), mixture = 1) %>% 
  set_engine("glmnet") %>% 
  set_mode("classification")

# Validação cruzada para ajuste do hiperparâmetro em 10 lotes
ridge_cv_split <- vfold_cv(treinamento_proc, v = 10)

doParallel::registerDoParallel() # Paralelizando os próximos comandos

# Definição do espaço de busca para o hiperparâmetro 'penalty'
val_gride_ridge <- grid_regular(penalty(range = c(0.001, 1)), levels = 20)

# Ajuste do hiperparâmetro usando validação cruzada
ridge_lambda_tune <- tune_grid(
  ridge,
  receita,
  resamples = ridge_cv_split,
  grid = val_gride_ridge,
  metrics = metric_set(accuracy),
  control = control_grid(save_pred = TRUE)
)

# Selecionando a melhor combinação de hiperparâmetros
ridge_best <- ridge_lambda_tune %>% 
  select_best("accuracy")

# Finalizando o modelo com os melhores hiperparâmetros
fit_ridge <- finalize_model(ridge, parameters = ridge_best) %>%
  fit(from_brazil ~ ., data = treinamento_proc)

# Previsões no conjunto de teste
ridge_fit_tunado <- predict(fit_ridge, teste_proc, type = "class") %>%
  bind_cols(teste_proc) %>%
  rename(Predito = .pred_class)

# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
ridge_fit_tunado <- ridge_fit_tunado %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))

9.2) Lasso Tunado

 # Definindo o modelo de regressão logística com penalidade Ridge
lasso <- logistic_reg(penalty = tune(), mixture = 0) %>% 
  set_engine("glmnet") %>% 
  set_mode("classification")

# Validação cruzada para ajuste do hiperparâmetro em 10 lotes
ridge_cv_split <- vfold_cv(treinamento_proc, v = 10)

doParallel::registerDoParallel() # Paralelizando os próximos comandos

# Definição do espaço de busca para o hiperparâmetro 'penalty'
val_gride_lasso <- grid_regular(penalty(range = c(0.001, 1)), levels = 20)

# Ajuste do hiperparâmetro usando validação cruzada
lasso_lambda_tune <- tune_grid(
  lasso,
  receita,
  resamples = ridge_cv_split,
  grid = val_gride_lasso,
  metrics = metric_set(accuracy),
  control = control_grid(save_pred = TRUE)
)

# Selecionando a melhor combinação de hiperparâmetros
lasso_best <- lasso_lambda_tune %>% 
  select_best("accuracy")

# Finalizando o modelo com os melhores hiperparâmetros
fit_lasso <- finalize_model(lasso, parameters = lasso_best) %>%
  fit(from_brazil ~ ., data = treinamento_proc)

# Previsões no conjunto de teste
lasso_fit_tunado <- predict(fit_lasso, teste_proc, type = "class") %>%
  bind_cols(teste_proc) %>%
  rename(Predito = .pred_class)

# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
lasso_fit_tunado <- lasso_fit_tunado %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))

9.3) Árvore de Decisào

# Definindo o modelo de árvore de decisão
arvore <- decision_tree(tree_depth = tune(), cost_complexity = tune()) %>% 
  set_engine("rpart") %>% 
  set_mode("classification")

# Paralelizando os próximos comandos
doParallel::registerDoParallel()

# Definição do espaço de busca para os hiperparâmetros 'tree_depth' e 'cost_complexity'
val_grid_arvore <- grid_regular(
  tree_depth(range = c(1, 10)), 
  cost_complexity(range = c(0.001, 0.1)), 
  levels = 10
)

# Criando o objeto de validação cruzada
cv_split <- vfold_cv(treinamento_proc, v = 10)

# Ajuste do hiperparâmetro usando validação cruzada
arvore_tune <- tune_grid(
  arvore,
  receita,
  resamples = cv_split,
  grid = val_grid_arvore,
  metrics = metric_set(accuracy)
)

# Coletando as métricas do modelo ajustado
arvore_metrics <- arvore_tune %>% 
  collect_metrics()

# Selecionando a melhor combinação de hiperparâmetros com base na acurácia
arvore_best <- arvore_tune %>% 
  select_best("accuracy")

# Finalizando o modelo com os melhores hiperparâmetros
fit_arvore_final <- finalize_model(arvore, parameters = arvore_best) %>%
  fit(from_brazil ~ ., data = treinamento_proc)

# Previsões no conjunto de teste
arvore_fit_tunado <- predict(fit_arvore_final, teste_proc, type = "class") %>%
  bind_cols(teste_proc) %>%
  rename(Predito = .pred_class)

# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
arvore_fit_tunado <- arvore_fit_tunado %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))

9.4) XBoost

# Define o modelo XGBoost com parâmetros para serem tunados
xgboost_model <- boost_tree(
  trees = tune(),           # Número de árvores
  tree_depth = tune(),      # Profundidade da árvore
  min_n = tune(),           # Número mínimo de observações
  loss_reduction = tune(),  # Redução de perda necessária para fazer uma divisão adicional
  mtry = tune(),            # Número de variáveis a considerar em cada divisão
  learn_rate = tune()       # Taxa de aprendizado
) %>%
  set_engine("xgboost") %>%
  set_mode("classification")


# Criando o objeto de validação cruzada
cv_split_boost <- vfold_cv(treinamento_proc, v = 10)


# Executa a sintonia (tuning) do modelo
xgboost_tune <- tune_grid(
  xgboost_model,
  receita,
  resamples = cv_split_boost,
  grid = 30,
  metrics = metric_set(accuracy)
)

## i Creating pre-processing data to finalize unknown parameter: mtry

# Coleta as métricas calculadas
xgboost_metrics <- xgboost_tune %>% 
  collect_metrics()

# Seleciona a melhor combinação de hiperparâmetros
xgboost_best <- xgboost_tune %>% 
  select_best("accuracy")

# Finaliza o modelo com os melhores hiperparâmetros
fit_xgboost <- finalize_model(xgboost_model, parameters = xgboost_best) %>%
  fit(from_brazil ~ ., data = treinamento_proc)

# Faz previsões no conjunto de teste
xgboost_fit_tunado <- predict(fit_xgboost, teste_proc, type = "class") %>% 
  bind_cols(teste_proc) %>%
  rename(Predito = .pred_class) %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))

9.5) Rede Neural

# Preparação dos dados para a Rede Neural **************************************

# Seleciona as variáveis preditoras e converte para matriz
X_trn <- treinamento_proc %>% 
  dplyr::select(-from_brazil) %>% 
  as.matrix()

X_tst <- teste_proc %>% 
  dplyr::select(-from_brazil) %>% 
  as.matrix()

# Normalização dos dados
X_trn <- scale(X_trn)
X_tst <- scale(X_tst, center = attr(X_trn, "scaled:center"), scale = attr(X_trn, "scaled:scale"))

# Converter as classes de destino para one-hot encoding
num_classes <- length(unique(treinamento_proc$from_brazil))
Y_trn <- to_categorical(as.numeric(factor(treinamento_proc$from_brazil)) - 1, num_classes)
Y_tst <- to_categorical(as.numeric(factor(teste_proc$from_brazil)) - 1, num_classes)

# Construção da Rede Neural ****************************************************

net <- keras_model_sequential() %>% 
  layer_dense(units = 64, activation = "relu", input_shape = ncol(X_trn)) %>% 
  layer_dense(units = 32, activation = "relu") %>% 
  layer_dense(units = 16, activation = "relu") %>% 
  layer_dense(units = num_classes, activation = "sigmoid")

net <- compile(
  net, 
  loss = "binary_crossentropy", 
  optimizer = "adam", 
  metrics = "accuracy"
)

summary(net)

## Model: "sequential"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  dense_3 (Dense)                    (None, 64)                      640         
##  dense_2 (Dense)                    (None, 32)                      2080        
##  dense_1 (Dense)                    (None, 16)                      528         
##  dense (Dense)                      (None, 2)                       34          
## ================================================================================
## Total params: 3282 (12.82 KB)
## Trainable params: 3282 (12.82 KB)
## Non-trainable params: 0 (0.00 Byte)
## ________________________________________________________________________________

# Treinamento da Rede Neural ****************************************************

history <- fit(net, X_trn, Y_trn, batch_size = 50, epochs = 30, validation_split = 0.2)

## Epoch 1/30
## 3/3 - 1s - loss: 0.7663 - accuracy: 0.4150 - val_loss: 0.8425 - val_accuracy: 0.4054 - 603ms/epoch - 201ms/step
## Epoch 2/30
## 3/3 - 0s - loss: 0.7310 - accuracy: 0.4898 - val_loss: 0.8241 - val_accuracy: 0.4324 - 49ms/epoch - 16ms/step
## Epoch 3/30
## 3/3 - 0s - loss: 0.7093 - accuracy: 0.5510 - val_loss: 0.8081 - val_accuracy: 0.5405 - 37ms/epoch - 12ms/step
## Epoch 4/30
## 3/3 - 0s - loss: 0.6980 - accuracy: 0.5646 - val_loss: 0.7919 - val_accuracy: 0.5405 - 36ms/epoch - 12ms/step
## Epoch 5/30
## 3/3 - 0s - loss: 0.6879 - accuracy: 0.5714 - val_loss: 0.7754 - val_accuracy: 0.5676 - 36ms/epoch - 12ms/step
## Epoch 6/30
## 3/3 - 0s - loss: 0.6790 - accuracy: 0.5918 - val_loss: 0.7614 - val_accuracy: 0.5135 - 36ms/epoch - 12ms/step
## Epoch 7/30
## 3/3 - 0s - loss: 0.6711 - accuracy: 0.6190 - val_loss: 0.7493 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 8/30
## 3/3 - 0s - loss: 0.6640 - accuracy: 0.6122 - val_loss: 0.7419 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 9/30
## 3/3 - 0s - loss: 0.6587 - accuracy: 0.6122 - val_loss: 0.7355 - val_accuracy: 0.4865 - 36ms/epoch - 12ms/step
## Epoch 10/30
## 3/3 - 0s - loss: 0.6532 - accuracy: 0.6259 - val_loss: 0.7338 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 11/30
## 3/3 - 0s - loss: 0.6486 - accuracy: 0.6395 - val_loss: 0.7335 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 12/30
## 3/3 - 0s - loss: 0.6446 - accuracy: 0.6395 - val_loss: 0.7306 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 13/30
## 3/3 - 0s - loss: 0.6416 - accuracy: 0.6395 - val_loss: 0.7322 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 14/30
## 3/3 - 0s - loss: 0.6392 - accuracy: 0.6531 - val_loss: 0.7338 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 15/30
## 3/3 - 0s - loss: 0.6362 - accuracy: 0.6599 - val_loss: 0.7352 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 16/30
## 3/3 - 0s - loss: 0.6344 - accuracy: 0.6667 - val_loss: 0.7365 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 17/30
## 3/3 - 0s - loss: 0.6321 - accuracy: 0.6599 - val_loss: 0.7387 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 18/30
## 3/3 - 0s - loss: 0.6306 - accuracy: 0.6667 - val_loss: 0.7434 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 19/30
## 3/3 - 0s - loss: 0.6282 - accuracy: 0.6599 - val_loss: 0.7497 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 20/30
## 3/3 - 0s - loss: 0.6259 - accuracy: 0.6667 - val_loss: 0.7515 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 21/30
## 3/3 - 0s - loss: 0.6242 - accuracy: 0.6667 - val_loss: 0.7589 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 22/30
## 3/3 - 0s - loss: 0.6232 - accuracy: 0.6735 - val_loss: 0.7651 - val_accuracy: 0.4595 - 34ms/epoch - 11ms/step
## Epoch 23/30
## 3/3 - 0s - loss: 0.6214 - accuracy: 0.6667 - val_loss: 0.7712 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 24/30
## 3/3 - 0s - loss: 0.6201 - accuracy: 0.6803 - val_loss: 0.7754 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 25/30
## 3/3 - 0s - loss: 0.6202 - accuracy: 0.6803 - val_loss: 0.7769 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 26/30
## 3/3 - 0s - loss: 0.6187 - accuracy: 0.6735 - val_loss: 0.7850 - val_accuracy: 0.4595 - 34ms/epoch - 11ms/step
## Epoch 27/30
## 3/3 - 0s - loss: 0.6173 - accuracy: 0.6735 - val_loss: 0.7917 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 28/30
## 3/3 - 0s - loss: 0.6166 - accuracy: 0.6735 - val_loss: 0.8002 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 29/30
## 3/3 - 0s - loss: 0.6177 - accuracy: 0.6735 - val_loss: 0.8054 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 30/30
## 3/3 - 0s - loss: 0.6165 - accuracy: 0.6735 - val_loss: 0.8137 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step

# Previsões e Avaliação *********************************************************

y_hat_net <- predict(net, X_tst)

## 3/3 - 0s - 58ms/epoch - 19ms/step

predicted_classes <- apply(y_hat_net, 1, which.max) - 1 # Obter as classes previstas

# Criando um dataframe temporário para as previsões da Rede Neural
previsoes_rede_neural <- cbind(teste_proc, Predito = as.factor(predicted_classes))

# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
previsoes_rede_neural <- previsoes_rede_neural %>%
  mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
         Predito = ifelse(Predito == 1, "Brasileiro", "Não Brasileiro"))

# 10) Matriz de Confusão

criar_df_confusao <- function(predicoes, observado, nome_modelo) {
  # Definindo os níveis na ordem desejada com "Brasileiro" primeiro
  niveis <- c("Brasileiro", "Não Brasileiro")
  
  # Convertendo para fatores com os níveis definidos
  predicoes <- factor(predicoes, levels = niveis)
  observado <- factor(observado, levels = niveis)
  
  cm <- confusionMatrix(predicoes, observado)
  df <- cbind(melt(as.matrix(cm$table)), Modelo = nome_modelo)
  return(df)
}

# Convertendo previsões e observações para fatores com os mesmos níveis
logistica$Predito <- factor(logistica$Predito)
logistica$Observado <- factor(logistica$Observado)
df_logistica <- criar_df_confusao(logistica$Predito, logistica$Observado, "Logística")

ridge_fit_tunado$Predito <- factor(ridge_fit_tunado$Predito)
ridge_fit_tunado$Observado <- factor(ridge_fit_tunado$Observado)
df_ridge <- criar_df_confusao(ridge_fit_tunado$Predito, ridge_fit_tunado$Observado, "Ridge Tunado")

lasso_fit_tunado$Predito <- factor(lasso_fit_tunado$Predito)
lasso_fit_tunado$Observado <- factor(lasso_fit_tunado$Observado)
df_lasso <- criar_df_confusao(lasso_fit_tunado$Predito, lasso_fit_tunado$Observado, "Lasso Tunado")

arvore_fit_tunado$Predito <- factor(arvore_fit_tunado$Predito)
arvore_fit_tunado$Observado <- factor(arvore_fit_tunado$Observado)
df_arvore <- criar_df_confusao(arvore_fit_tunado$Predito, arvore_fit_tunado$Observado, "Árvore de Decisão")

xgboost_fit_tunado$Predito <- factor(xgboost_fit_tunado$Predito)
xgboost_fit_tunado$Observado <- factor(xgboost_fit_tunado$Observado)
df_xgboost <- criar_df_confusao(xgboost_fit_tunado$Predito, xgboost_fit_tunado$Observado, "XGBoost")

# Para a Rede Neural, ajustando as previsões
predicted_classes_net <- ifelse(previsoes_rede_neural$Predito == 1, "Brasileiro", "Não Brasileiro")
predicted_classes_net <- factor(predicted_classes_net)
previsoes_rede_neural$Observado <- factor(previsoes_rede_neural$Observado)
df_rede_neural <- criar_df_confusao(predicted_classes_net, previsoes_rede_neural$Observado, "Rede Neural")

# Combinando todos os dataframes
df_all <- rbind(df_logistica, df_ridge, df_lasso, df_arvore, df_xgboost, df_rede_neural)


ggplot(data = df_all, aes(x = Prediction, y = Reference)) +
  geom_tile(aes(fill = value), color = 'darkblue') +
  geom_text(aes(label = sprintf("%d", value)), vjust = 1, size = 3) +
  scale_fill_gradient(low = "lightblue", high = "steelblue") +
  scale_y_discrete(limits = rev(levels(df_all$Reference))) +  # Isso irá inverter a ordem no eixo Y
  theme_minimal() +
  labs(title = "Comparação de Matrizes de Confusão", x = "Predição", y = "Real", fill = "Contagem") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text()) +
  facet_wrap(~ Modelo, ncol = 3)

10) Curva RoC

# Para o modelo Logistica
logistica_probs <- predict(logistica_simples, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo Ridge
ridge_probs <- predict(fit_ridge, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo Lasso
lasso_probs <- predict(fit_lasso, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo de Árvore de Decisão
arvore_probs <- predict(fit_arvore_final, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo XGBoost
xgboost_probs <- predict(fit_xgboost, teste_proc, type = "prob")[[".pred_1"]]
# Para a Rede Neural
rede_neural_probs <- predict(net, X_tst)[, 1]

## 3/3 - 0s - 12ms/epoch - 4ms/step

resultados_preds <- bind_rows(
  tibble(
    modelo = "Logística", 
    marcador = logistica_probs, 
    resposta = as.factor(teste$from_brazil)
  ), 
  tibble(
    modelo = "Ridge", 
    marcador = ridge_probs,
    resposta = as.factor(teste$from_brazil)
  ), 
  tibble(
    modelo = "Lasso",
    marcador = lasso_probs, 
    resposta = as.factor(teste$from_brazil)
  ), 
  tibble(
    modelo = "Árvore de Decisão",
    marcador = arvore_probs,
    resposta = as.factor(teste$from_brazil)
  ),
  tibble(
    modelo = "XGBoost",
    marcador = xgboost_probs,
    resposta = as.factor(teste$from_brazil)
  ),
  tibble(
    modelo = "Neural",
    marcador = rede_neural_probs,
    resposta = as.factor(teste$from_brazil)
  )
)

# Calculando a curva ROC para cada modelo e plotando
roc_plot <- resultados_preds %>% 
  group_by(modelo) %>%
  roc_curve(resposta, marcador, event_level = "second") %>% 
  autoplot() +
  theme(legend.position = "right")

print(roc_plot)

11) Analisando as Métricas

( Acurácia, Precisão , Recall e F1 - Score)

#Avaliando metricas _ 

lista_modelos <- list(df_logistica, df_ridge, df_lasso, df_arvore, df_xgboost, df_rede_neural)

# Calculando as métricas diretamente dos dados fornecidos
metricas_modelos <- lapply(lista_modelos, function(df) {
  verdadeiro_positivo <- df$value[df$Prediction == "Brasileiro" & df$Reference == "Brasileiro"]
  falso_positivo <- df$value[df$Prediction == "Não Brasileiro" & df$Reference == "Brasileiro"]
  falso_negativo <- df$value[df$Prediction == "Brasileiro" & df$Reference == "Não Brasileiro"]
  verdadeiro_negativo <- df$value[df$Prediction == "Não Brasileiro" & df$Reference == "Não Brasileiro"]
  
  acuracia <- (verdadeiro_positivo + verdadeiro_negativo) / sum(df$value)
  precisao <- verdadeiro_positivo / (verdadeiro_positivo + falso_positivo)
  recall <- verdadeiro_positivo / (verdadeiro_positivo + falso_negativo)
  f1 <- ifelse((precisao + recall) == 0, 0, 2 * (precisao * recall) / (precisao + recall))
  
  tibble(
    Modelo = unique(df$Modelo),
    Acuracia = acuracia,
    Precisao = precisao,
    Recall = recall,
    F1 = f1
  )
}) %>% bind_rows()

# Tratando NA, convertendo métricas para porcentagem e ordenando por acurácia
metricas_modelos <- metricas_modelos %>%
  mutate(
    across(c(Acuracia, Precisao, Recall, F1), ~ifelse(is.na(.), 0, round(. * 100, 0)))
  ) %>%
  arrange(desc(Acuracia))

# Formatação dos dados para a visualização com gt
gt_table <- metricas_modelos %>%
  gt() %>%
  tab_header(
    title = "Métricas de Desempenho dos Modelos"
  ) %>%
  cols_label(
    Modelo = "Modelo",
    Acuracia = "Acurácia",
    Precisao = "Precisão",
    Recall = "Recall",
    F1 = "F1 Score"
  ) %>%
  fmt_number(
    columns = c("Acuracia", "Precisao", "Recall", "F1"),
    decimals = 0,
    pattern = "{x}%"
  ) %>%
  tab_options(
    heading.background.color = "gray",
    column_labels.font.size = 12,
    data_row.padding = px(10)
  ) %>%
  tab_style(
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      columns = c("Acuracia", "Precisao", "Recall", "F1")
    )
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels(
      columns = c("Modelo", "Acuracia", "Precisao", "Recall", "F1")
    )
  )

# Exibindo a tabela
gt_table

Modelo	Acurácia	Precisão	Recall	F1 Score
Métricas de Desempenho dos Modelos
XGBoost	65%	65%	62%	63%
Rede Neural	54%	0%	0%	0%
Logística	50%	62%	47%	53%
Lasso Tunado	48%	81%	46%	59%
Ridge Tunado	46%	100%	46%	63%
Árvore de Decisão	46%	100%	46%	63%

12) Analise Não Supervisionada - PCA

A ánalise foi feita para os países mais avaliados em número de observacões dentro das observacões ( Top 10)

12.1) Preparando PCA

# Passando os dados tratados para o nova variavel
dados_clust <- dados_final

# Contando os registros para cada país para levar somente o top 10 mais avaliados em frequencia
country_counts <- dados_clust %>%
  count(country_of_origin) %>%
  arrange(desc(n))

# Selecionando os top 6 países com mais registros
top_countries <- head(country_counts, 10)

# Filtrando 'dados_clust' para incluir apenas os top 6 países
filtered_clust <- dados_clust %>%
  filter(country_of_origin %in% top_countries$country_of_origin)

# Agrupando e somando os atributos para os países filtrados
df <- filtered_clust %>%
  group_by(country_of_origin) %>%
  summarise(across(c(aroma, flavor, aftertaste, acidity, body, balance, uniformity, clean_cup, sweetness), ~sum(., na.rm = TRUE)))

# Definir a coluna 'regiao' como row names e remover a coluna do dataframe

df <- df %>%
  column_to_rownames(var = "country_of_origin")


X <- scale(df, 
           center = TRUE, # centraliza os dados
           scale = TRUE) # escalona os dados (pois estao em medidas diferentes)

12.2) Aplicando PCA e Analisando

pca <- prcomp(X) # aplica o PCA

# Todos os atributos para PC1 tem o mesmo peso nao mostrando atributos que se sobressaem

pca$rotation <- -pca$rotation # troca o sinal das cargas
pca$x <- -pca$x # troca o sinal dos scores

Phi <- pca$rotation # matriz de cargas
head(Phi)

##                  PC1        PC2         PC3         PC4        PC5         PC6
## aroma      0.3333679 0.01157345  0.08366266 -0.04978241  0.4417232 -0.15801712
## flavor     0.3333591 0.16693211  0.16272529 -0.05493036 -0.6843610  0.14748734
## aftertaste 0.3333299 0.34140572 -0.37119469 -0.36685208 -0.1632435 -0.26629708
## acidity    0.3333312 0.01458445  0.79592318  0.14359109  0.1305645  0.01362871
## body       0.3333602 0.18745743  0.10937622 -0.12224716 -0.1729299 -0.06988413
## balance    0.3333115 0.50462916 -0.23607065  0.12983608  0.4671684  0.30713836
##                    PC7         PC8         PC9
## aroma      -0.17416192  0.61170571  0.50462770
## flavor     -0.08661671  0.51395144 -0.26411654
## aftertaste -0.47601668 -0.38692467  0.16156189
## acidity    -0.26335184 -0.37533772 -0.08426192
## body        0.77360425 -0.20443827  0.38958452
## balance     0.16453901  0.02848641 -0.47032037

Z <- pca$x # matriz de scores
head(Z)

##                  PC1          PC2          PC3          PC4           PC5
## Brazil      1.096309  0.007548140 -0.011464123 -0.004428219 -0.0152341163
## Colombia    3.347607  0.057000370 -0.031358419  0.007881091  0.0054410453
## Costa Rica -2.357759 -0.003584463 -0.003792201 -0.003037407  0.0083292517
## Ethiopia   -2.594418  0.031025358  0.019949117 -0.006628561 -0.0043695017
## Guatemala   3.149011 -0.006995421  0.027747603  0.038912421  0.0005027964
## Honduras   -2.365499 -0.033384541 -0.003929353  0.006117134 -0.0008986453
##                     PC6           PC7           PC8           PC9
## Brazil      0.002174996  0.0001237894 -0.0005462686 -2.258589e-05
## Colombia    0.001064819  0.0002477694  0.0001204255  1.233801e-05
## Costa Rica  0.001235835 -0.0035754693  0.0012226895 -3.276435e-05
## Ethiopia    0.001609311 -0.0046573502  0.0003097058  2.670071e-05
## Guatemala  -0.000535794 -0.0001326772 -0.0001578867 -5.862193e-06
## Honduras    0.003455312  0.0049293446  0.0018220648  1.308488e-05

#Analises PCA_______________-

# Resumo do PCA para verificar a variação explicada por cada componente principal
summary(pca)

## Importance of components:
##                           PC1     PC2     PC3     PC4      PC5      PC6
## Standard deviation     2.9997 0.03520 0.01884 0.01808 0.006653 0.003218
## Proportion of Variance 0.9998 0.00014 0.00004 0.00004 0.000000 0.000000
## Cumulative Proportion  0.9998 0.99992 0.99996 0.99999 1.000000 1.000000
##                             PC7     PC8       PC9
## Standard deviation     0.003044 0.00126 1.786e-05
## Proportion of Variance 0.000000 0.00000 0.000e+00
## Cumulative Proportion  1.000000 1.00000 1.000e+00

# Obtendo as cargas para cada atributo nos componentes principais
loadings <- abs(pca$rotation) 

# Para identificar os atributos mais importantes para o primeiro componente principal (PC1)
atributos_importantes_pc1 <- sort(loadings[,1], decreasing = TRUE)

# Imprimindo os atributos mais importantes para PC1
print(atributos_importantes_pc1)

##      aroma       body     flavor  clean_cup    acidity uniformity aftertaste 
##  0.3333679  0.3333602  0.3333591  0.3333313  0.3333312  0.3333308  0.3333299 
##    balance  sweetness 
##  0.3333115  0.3332782

# Convertendo a matriz Phi em uma tibble, preservando os nomes das linhas
contribuicoes <- as_tibble(Phi, rownames = "country_of_origin") %>%
  mutate(Contribuicao_PC1 = (PC1^2) * 100 / sum(PC1^2)) %>%
  dplyr::select(country_of_origin, Contribuicao_PC1)

# Visualizando as contribuições
print(contribuicoes)

## # A tibble: 9 × 2
##   country_of_origin Contribuicao_PC1
##   <chr>                        <dbl>
## 1 aroma                         11.1
## 2 flavor                        11.1
## 3 aftertaste                    11.1
## 4 acidity                       11.1
## 5 body                          11.1
## 6 balance                       11.1
## 7 uniformity                    11.1
## 8 clean_cup                     11.1
## 9 sweetness                     11.1

12.3) Análise Gráfica

# o grafico abaixo mostra o percentual explicado da variancia de cada componente
fviz_eig(pca, addlabels = TRUE) + 
  labs(x = "Componente Principal",
       y = "Percentual explicado da variância")

# Grafico mostra a PCA explicando quase toda a variancia da base

# biplot do factoextra
fviz_pca_biplot(pca, 
                repel = TRUE, 
                xlab = "PC1 - Intensidade do Aroma", 
                ylab = "PC2 - Qualidade do Sabor",
                labelsize = 3, 
                geom = c("point", "text"), 
                addEllipses = FALSE) + 
  theme(plot.margin = margin(0, 0, 0, 0, "cm"), 
        text = element_text(size = 10)) +
  theme(axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12), 
        axis.text.x = element_text(size = 8),  
        axis.text.y = element_text(size = 8))

12.4) Análise Brasil PCA

Brasil: Está próximo do centro no eixo PC1, indicando uma intensidade de aroma moderada. Os atributos “body”, “balance” e “aftertaste” apontam na direção do Brasil, sugerindo que esses atributos são fortes no café brasileiro.

13) Conclusão

Nossa expectativa era alcançar uma precisão superior a 65% na identificação do café brasileiro. Contudo, enfrentamos desafios significativos devido ao desbalanceamento dos dados. Esse obstáculo nos proporcionou um valioso aprendizado sobre como manejar eficientemente tais problemas. Empregamos e descartamos múltiplos modelos de regressão altamente colinearizados, e adotamos diversas abordagens e técnicas para mitigar essas dificuldades. Essa jornada nos levou à conclusão de que uma acurácia de 65% e a utilização do método XGBoost representam a melhor performance possível para determinar se um café é brasileiro ou não.

PADS - Aprendizado de Máquina 2Projeto Final - Café Brasileiro

Bruno Mernick, Giovanna Rodrigues e Mariana Aukar

02/12/2023